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Lattice calculation of non-Gaussianity from preheating 
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If light scalar fields are present at the end of inflation, their non-equilibrium dynamics such as 
parametric resonance or a phase transition can produce non-Gaussian density perturbations. We 
show how these perturbations can be calculated using non-linear lattice field theory simulations and 
the separate universe approximation. In the massless preheating model, we find that some parameter 
values are excluded while others lead to acceptable but observable levels of non-Gaussianity. This 
shows that preheating can be an important factor in assessing the viability of inflationary models. 
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PACS numbers: 98.80.Cq, 11.15.Kc 

Predicting the observed nearly scale-invariant and 
Gaussian primordial density perturbations has been the 
greatest success of the slow-roll inflationary paradigm. 
However, these predictions are very generic and shared 
by a vast number of different inflationary models. Even 
the small observed deviation from scale invariance can be 
expressed purely in terms of slow roll parameters, and is 
therefore not enough to distinguish between models. 

In contrast, if deviation from Gaussianity is observed, 
it will at least in principle provide much more information 
and could allow us to single out the correct model. The 
level of non-Gaussianity is commonly characterised using 
the parameter /nl, which was originally defined [l[ for a 
specific "local" type of non-Gaussianity by 
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where C, is the curvature perturbation and Co is a Gaus- 
sian random field. The definition of /nl has since been 
extended to cover arbitrary non-Gaussian fields 0,0], in 
which case it also becomes scale-dependent. Even this 
generalised /nl does not give a full description of the 
non-Gaussianity, but it is nevertheless a useful quantity. 

Observational constraints on /nl are currently not par- 
ticularly strong |/nl| ^ 100, but they are Hkely to 
improve significantly with the Planck satellite. Unfortu- 
nately, the level of non-Gaussianity produced during in- 
flation in single-field models is very small 0], |/nl| ^ 1, 
and is unlikely to be ever observed. 

In many inflationary models the end of inflation in- 
volves non-linear non-equilibrium dynamics, such as a 
phase transition or parametric resonance. It is natural to 
ask whether these processes can produce density pertur- 
bations and how non-Gaussian this contribution would 
be (see recent review [5] and references therein). It is 
clear that causality prohibits generation of super-Hubble 
perturbations from nothing. However, if there is a light 
scalar field, it would acquire similar scale invariant per- 
turbations during inflation to the inflaton field, and these 



may well be converted into density perturbations at the 
end of inflation. 

In this letter, we show how the perturbations produced 
by non-equilibrium dynamics at the end of inflation can 
be calculated in a non-linear way using lattice field the- 
ory simulations. The method is very general and can be 
applied to any inflationary model in which the dynamics 
are dominated by bosonic fields. For concreteness, we 
will focus on a particular inflationary model known as 
masslesspreheating [1, 0, Q • Recent approximate calcu- 
lations [3, [lO] have suggested that preheating could pro- 
duce very large levels of non-Gaussianity in this model, 
and we show that this is indeed the case for large parts 
of the parameter space. We will also discuss at a general 
level the conditions required for effective generation of 
perturbations at the end of inflation. 

The model consists of an inflaton field (j) coupled to a 
massless scalar field x? with the potential 
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During inflation, x is approximately zero, and the model 
behaves the same way as the ordinary single-field cj)'^ 
chaotic inflation model. If g^/X < 0(1), the mass 
™x ~ 9't' of A field is less than the Hubble rate, 
and therefore the x field acquires a similar scale-invariant 
spectrum of perturbations to the inflaton field (j). 

When inflation ends, the inflaton starts to oscil- 
late with a decreasing amplitude. In terms of rescaled 
field (j) = a(j) and rescaled conformal time r defined by 
dr = a~^X^^'^4>inidt, the oscillations are approximately 
described by the Jacobi cosine function [11], 0(t) = 
(^inicn(r, l/\/2). This gives rise to an oscillatory mass 
term for the x field. At linear level, a Fourier mode of 
the rescaled field x — ^x with wave number k satisfies 
the Lame equation. 
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In the space of the two constant parameters and 
g^/X, Eq. ^ has resonance bands in which the solu- 



tion grows exponentially, Xk{T) 



-"■^/(r), where fi is 



2 



known as the Floquet index and /(r) is a periodic func- 
tion (see Fig. 4 in Ref. Q)- This means that energy is 
transferred rapidly from the inflaton (j) to the x held, 
which is known as preheating Ts*!. The resonance 
bands stretch through the (g^/A, k^) plane diagonally so 
that for any value of /X there are some resonant modes, 
but when <C A, the bands are narrow and the res- 
onance is weak. Effective preheating therefore requires 
.9VA > 1. 

This model has been studied extensively using lattice 
field theory simulations [1, , but the focus has always 
been on field dynamics rather than on its effect on cur- 
vature perturbations. In a recent work , gravitational 
effects were introduced at linearised level by coupling the 
fields to metric perturbations. This way, one can study 
their back- reaction to the field dynamics but, because 
only linear small-scale perturbations are included, not 
observable length scales or non-Gaussianity. 

Our approach combines lattice field theory with the 
widely used separate universe approximation [l^. It 
states that points in space separated by more than a 
Hubble distance cannot interact and will therefore evolve 
independently of each other. As long as each Hub- 
ble volume is approximately isotropic and homogeneous, 
one can approximate them by separate Friedmann- 
Robertson- Walker (FRW) universes. When the universe 
has reached equilibrium, the curvature perturbation C 
is equal to the perturbation of the logarithm of the 
scale factor a on a constant energy density hypersurface, 
C = 6N = 6\iia\H, which one can find by solving the 
Friedmann equation for each universe separately. 

At the end of inflation, all fields except light scalar 
fields vanish. It is therefore sufficient to consider the de- 
pendence of N on the "initial" values of the light scalars 
(pi and Xi by which we mean their average values in the 
given Hubble volume at the end of inflation. The model is 
symmetric under Xi ~^ ~X'n which implies dN/dxi = 0. 
We can therefore Taylor expand 
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where we have subtracted the averages {S(t)f) and {Sxf) 
to keep the unperturbed value of C zero. Assuming, as 
we do, that the initial field values are Gaussian, the first 
term described the usual Gaussian contribution, and the 
second term is a local non-Gaussian contribution that is 
present already in single-field models. The third term, 
however, is more interesting, since it represents the con- 
tribution from the other light field x- To measure the 
level of this contribution, we need to find the dependence 
of = In a on the initial value of x- 

Previous applications of the separate universe ap- 
proach they have generally assumed that not only the 
metric but also the field confi gura tions are homogeneous 
inside each separate universe (l7l . [l8| . the dynamics are 



then described by a system of coupled ODEs, 
• dV 

■ dV 

X + 3Hc^+— = 0, 
dx 



(5) 



where p = ^cjp' -\- \x^ -\- V{(j), x) is the energy density. 
Solving these equations and calculating the scale factor 
a at some suitably chosen H gives a non-linear way of 
calculating the curvature perturbation C,. 

The contribution by the x held to the non-Gaussianity 
is not of the simple "local" type ^ , but defining /nl as a 
suitable ratio of the three-point and two-point correlation 
functions of C, Boubekeur and Lyth obtained for this type 
of non-Gaussianity an effective value of Q 
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where we have used the expressions = (i?^/47r^) and 
Vq — (V"/247r^eMp[) for the power spectra of x and 
The logarithm reflects an infrared divergence, which is 
cut off by the length scale at which the averages in (6) 
are computed, i.e., the maximum observable scale. This 
expression is only valid if the non-Gaussianity is small, 
/nl ^ 10^. In any case, the observational limits on /nl 
lead to the constraint 
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The key finding in earlier work on this model [ITI . |l8| 
was that the dependence of N on the initial x appears 
random, suggesting chaotic dynamics. This would mean 
that the the contribution to the curvature perturbation 
is simply white noise, with no observable effects on large 
scales. As we will see, this chaotic behaviour is, in 
fact, merely a consequence of ignoring the inhomogeneous 
modes completely. When they are included, individual 
degrees of freedom may still behave chaotically, but the 
averaged dynamics, which determines the expansion of 
the universe and thereby the curvature perturbation, is 
non-chaotic. 

To include the effects of sub-Hubble inhomogeneities, 
we employ a hierarchy between the two relevant length 
scales. The length scales relevant for the non-equilibrium 
dynamics, which we call microscopic, are smaller than 
the horizon size, but we are interested in the curvature 
perturbation on macroscopic, super-Hubble scales that 
are observable today. 

We approximate the microscopic dynamics by defin- 
ing the inflaton and the other relevant fields as classical 
fields on a comoving lattice. The lattice size is chosen 
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to be large enough to describe the microscopic dynamics 
accurately but smaller than the horizon at all times, so 
that the space-time geometry inside the lattice can be 
well approximated by the FRW metric. The simulations 
in Ref. 15| confirm that sub-horizon metric perturba- 
tions can be safely ignored. We can therefore use the 
standard classical field theory methods that have been 
used to study preheating over the last decade. 

Since the lattice size is larger than the characteristic 
length scale of the microscopic dynamics, the average en- 
ergy density is the whole Hubble volume is well approx- 
imated by the average energy density p in our lattice, 
given by 
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where the integration is over the lattice and L is the co- 
moving lattice size. 

We can now use p in the Friedmann equation to cal- 
culate the evolution of the scale factor. In a sense, we 
are therefore using the field theory simulation to calcu- 
late the equation of state inside the Hubble volume. This 
same approach has been used in the past to study the for- 
mation of primordial black holes during preheating [l9(. 
The field evolution is given by the full classical equation 
of motion. 
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We solved this system of equations in conformal time 
{drj = a~^dt) using the second-order Runge-Kutta algo- 
rithm for the field equations coupled to an Euler method 
for the Friedmann equation 20]. 

To set up the initial conditions, we followed the stan- 
dard approach in studies of preheating [1, [2l| . We drew 
the initial values of the fields and their time derivatives 
from a Gaussian ensemble that has the same two-point 
functions as the quantum mechanical vacuum state. 



iXkXq) 
iXkXq) 



(2^)35(fc + g) 
{2TTfS{k + q) 



1 

2ujk ' 
2 ' 



(10) 



with all other two-point correlators vanishing and the 
(/i-field homogeneous at the beginning of the simulation. 
In the linearised classical theory, these two-point func- 
tions evolve exactly the same way as in the quantum 
theory. On the other hand, when the dynamics become 
non-linear, occupation numbers of the modes are gener- 
ally large, and the classical field equations should again 
describe the dynamics very well. For each set of param- 
eters, we repeated the simulation between 60 and 240 




FIG. 1: The dependence of A'^ = In a on Xi measured at H ■ 
5.53X 10"^^Mpi for g'^/X = 1.192 (unfilled circles) and g^/X ■ 
1.875 (filled circles). The curves are quadratic fits (|lip . 



times with different random realisation of the initial fluc- 
tuations. 

We fixed the infiaton self-coupling to A = 7 x 10^^^, 
and started the simulation at = 5 in Planck units. We 
used lattices of up to 32'^ points with spacing 6x — 1.25 x 
10^ and time step St = 4 x 10^. For these parameters, 
the contribution to the energy density by the "quantum" 
fluctuations (fTU|) . pfluct ~ Sx~'^ ~ 10~^^, is well below 
the physical energy density Pphys ~ ^'t'fni ^ 10~^^, and 
should not affect the dynamics. 

To determine the curvature perturbation C we have to 
measure the dependence of the number of e-foldings N = 
In a on the initial value of x at a suitable fixed value of 
the Hubble rate such that the system has reached a quasi- 
equilibrium state. Figure [T] shows A^(x) measured at H = 
5.53X lO-i^Mpi for g^/X = 1.192 and g^/X = 1.875. The 
data appear smooth, within errors, in contrast with the 
previous separate universe calculations [l3, [3 . There is 
significant scatter between individual runs, but since that 
variation is random and uncorrelated, it will contribute 
a white noise signal that is unobservable on large scales. 

According to Eq. ([4|) , the new contribution to the cur- 
vature perturbation depends on the second derivative 
d^N/dXi where Xi is the value of the zero mode at the 
beginning of the simulation. We measure it by fitting the 
data at low Xi with a quadratic function. 
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so that d^N/dxf = 2c. The curves in Fig. [T] show the 
fits for g'^/X = 1.192 and g^/X = 1.875. 

Table U shows the fitted values of d'^N/dxf for four 
different couplings g^/X. These values do not depend 
significantly on the value of H at which they are mea- 
sured. At g^/X — 1.05 the measured value corresponds 
to /nl ~ O (1). This is within observational bounds ([7]) 
and would be observable in future experiments, whereas 
the three other points are excluded by current observa- 
tions. In fact, they are so large that the contribution from 
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1.192 
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0.157 


1.875 


■|^gll.26±0.05 


0.237 


2.7 
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TABLE I: The second derivative d^N/dXi determined from 
the fit ([11| at ii" = 5.53 x 10"^^Mpi for different values of 
g^/\- The third column shows the Floquet index of the ho- 
mogeneous mode. 

preheating dominates the curvature perturbation. Even 
the spectrum would therefore disagree with observations. 

Comparison of the three largest values of g^/A in Ta- 
bic |T] demonstrates the importance of the inhomogeneous 
modes for the effect. The resonance is at its strongest at 
.g^/A = 1.875, and earlier studies [13, EE El which ig- 
nored inhomogeneous modes would have suggested that 
effect would be strongest at this parameter value. Like- 
wise, they would have suggested that the parameter val- 
ues g^/X = 2.7 and g^/X = 1.192 would have produced 
a smaller effect because the resonance is weaker. How- 
ever, our results show strongest effect at g^/X — 2.7 and 
weaker at g^/X = 1.192. 

This finding has a simple explanation. What matters 
is not the absolute strength of the resonance, but the 
relative growth rates of the zero mode and the inhomo- 
geneous modes. Near the lower edge of the resonance 
band g^/A = I, the growth rate of the zero mode is very 
small but inhomogeneous modes grow with a significant 
rate. Hence, d^N/dxf is very small. As g^/X is increased, 
the growth rate of the zero mode grows, and d'^N/dxf 
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grows very rapidly. At 5^/A « 1.4, the two growth rates 
are equal, and the dependence on g^/X is weaker. Above 
this value, d'^N/dxf starts to grow more rapidly again, 
as the growth rate of the inhomogeneous modes starts to 
fall. This agrees with our findings, and it will be very 
interesting to map this behaviour out including the sign 
of d'^N/dxf in detail in future studies ^2^. 

To summarise, we have demonstrated a method for 
computing the curvature perturbation non-linearly. The 
method can be readily applied to any inflationary model 
in which the dynamics are dominated by bosonic fields, 
ft requires an average over random initial conditions and 
therefore the results will inevitably have statistical errors, 
but as our results show, it is possible to achieve reason- 
ably good signal-to-noise ratios with existing computing 
technology. 

The strength of the mechanism that generates the cur- 
vature perturbations during preheating depends on how 
sensitive the expansion of the universe is on the initial 
value of the homogeneous mode. This is determined by 
the relative growth rates of the zero mode and the in- 
homogeneous modes. Thus, calculations that ignore the 
effect of inhomogeneous modes cannot be trusted. 

Our results rule out parts of the parameter space in 
the massless preheating model, because they would lead 
to extremely high level of non-Gaussian curvature per- 
turbations, but they also show that others are compati- 
ble with observations. They demonstrate that reheating 
mechanisms can be crucial for the viability of inflationary 
models. 
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Our simulation code, which was used to produce the numerical results in Fig. 1 and Table 1, contained an error. 
The lattice discretization of the expression (8) for the energy density was incorrect, and the calculated quantity was 
therefore not the conserved energy that corresponds to the discretized lattice field equations. As a consequence the 
numerical results cannot be trusted. This does not affect the validity of the method we described in our paper. 
Corrected simulation results will appear in Ref. [l| . We thank the authors of that paper for pointing out this error to 
us. 
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